This notebook re-implements the code from D. Higham's paper Modeling and Simulating Chemical Reactions.
This is a simple implementation of the Euler-Maruyama method to simulate the Chemical Langevin Equation for the Michelis-Menten system.
In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import ode
In [2]:
# Stoichiometric matrix
V = np.array([[-1.0, 1.0, 0.0],[-1.0, 1.0, 1.0],[1.0, -1.0, -1.0],[0.0, 0.0, 1.0]])
# Parameters and Initial Conditions
nA = 6.023e23 # Avagadro's number
vol = 1e-15 # volume of system
Y = np.zeros((4,))
c = np.zeros((3,))
d = np.zeros_like(c)
a = np.zeros_like(c)
Y[0] = round(5e-7 * nA * vol) # molecules of substrate
Y[1] = round(2e-7 * nA * vol) # molecules of enzyme
c[0] = 1.0e6 / (nA * vol)
c[1] = 1.0e-4
c[2] = 0.1
tfinal = 50.0
L = 250
tau = tfinal / L # stepsize
Yvals = np.zeros((4,L+1))
Yvals[:, 0] = Y
for k in range(L):
a[0] = c[0] * Y[0] * Y[1]
a[1] = c[1] * Y[2]
a[2] = c[2] * Y[2]
for i in range(3):
d[i] = tau * a[i] + np.sqrt(abs(tau * a[i])) * np.random.randn()
Y = Y + d[0] * V[:, 0] + d[1] * V[:, 1] + d[2] * V[:, 2]
Yvals[:, k+1] = Y
tvals = np.linspace(0.0, tfinal, L+1)
plt.plot(tvals, Yvals[0,:], 'go-', label="Substrate")
plt.plot(tvals, Yvals[3, :], 'r*-', label="Products")
plt.xlim((0.0, 55.0))
plt.ylim((0.0, 310.0))
plt.xlabel("Time", size=14)
plt.ylabel("Molecules", size=14)
plt.legend(fontsize=14)
plt.show();
Still simulating the same problem, but computing the reaction rates.
In [3]:
def mm_rre(t, y, k):
yprime = np.zeros_like(y)
yprime[0] = -k[0] * y[0] * y[1] + k[1] * y[2]
yprime[1] = -k[0] * y[0] * y[1] + (k[1] + k[2]) * y[2]
yprime[2] = k[0] * y[0] * y[1] - (k[1] + k[2]) * y[2]
yprime[3] = k[2] * y[2]
return yprime
In [4]:
tspan = np.array([0.0, 50.0])
yzero = np.array([5.0e-7, 2.0e-7, 0.0, 0.0])
k = np.array([1.0e6, 1.0e-4, 0.1])
r = ode(mm_rre).set_integrator('vode', method='bdf', order=15)
r.set_f_params(k)
r.set_initial_value(yzero, tspan[0])
dt = 1.0
y = np.zeros((4,int(tspan[1]/dt)+1))
y[:,0] = yzero
k = 1
while r.successful() and r.t < tspan[1]:
r.integrate(r.t + dt)
y[:,k] = r.y
k+=1
In [5]:
tvals = np.linspace(tspan[0], tspan[1], int(tspan[1]/dt)+1)
plt.plot(tvals, y[0, :], 'go-', markersize=10, label='Substrate')
plt.plot(tvals, y[3, :], 'r*-', markersize=10, label='Product')
plt.xlabel("Time", size=14)
plt.ylabel("Concentration", size=14)
plt.legend(fontsize=12)
plt.show()
A simple implementation of Gillespie's algorithm, or the Stochastic Simulation Algorithm, applied to the same system as above.
In [6]:
# Stoichiometric matrix
V = np.array([[-1.0, 1.0, 0.0],[-1.0, 1.0, 1.0],[1.0, -1.0, -1.0],[0.0, 0.0, 1.0]])
# Parameters and Initial Conditions
nA = 6.023e23 # Avagadro's number
vol = 1e-15 # volume of system
X = np.zeros((4,))
c = np.zeros((3,))
d = np.zeros_like(c)
a = np.zeros_like(c)
X[0] = round(5e-7 * nA * vol) # molecules of substrate
X[1] = round(2e-7 * nA * vol) # molecules of enzyme
c[0] = 1.0e6 / (nA * vol)
c[1] = 1.0e-4
c[2] = 0.1
t = 0.0
tfinal = 50.0
count = 1
tvals = [0.0]
Xvals = [list(X)]
while t < tfinal:
a[0] = c[0] * X[0] * X[1]
a[1] = c[1] * X[2]
a[2] = c[2] * X[2]
asum = np.sum(a)
j = np.argmax(np.random.rand() < np.cumsum(a / asum))
tau = np.log(1.0 / np.random.rand()) / asum
X += V[:, j]
t += tau
count += 1
tvals.append(t)
Xvals.append(list(X))
In [7]:
L = len(tvals)
tnew = np.zeros(2*L-1)
tnew[1:-1:2] = tvals[1:]
tnew[2::2] = tvals[1:]
tnew[0] = tvals[0]
Svals = np.array(Xvals)[:, 0]
ynew = np.zeros(2*L-1)
ynew[::2] = Svals
ynew[1:-1:2] = Svals[:-1]
plt.plot(tnew, ynew, 'go-', label = "Substrate")
Pvals = np.array(Xvals)[:, 3]
ynew = np.zeros(2*L-1)
ynew[::2] = Pvals
ynew[1:-1:2] = Pvals[:-1]
plt.plot(tnew, ynew, 'r*-', label = "Product")
plt.xlabel("Time", size=14)
plt.ylabel("Molecules", size=14)
plt.xlim((0.0, 55.0))
plt.ylim((0.0, 310.0))
plt.legend(fontsize=14)
plt.show()